

simCohort = function(parameters, Ncomm=Ngroups[Dgroup], varying) {                   



# get random effect
for (Dvarying in names(varying))
    parameters[[Dvarying]] = varying[[Dvarying]][1]

NindivPerCommunity = round(parameters$size / Ncomm)
parameters$size = Ncomm * NindivPerCommunity


result = data.frame(
  group = rep(1:Ncomm, NindivPerCommunity)
)


# csd level covariates, eg, law enforcement


groupsWithEffect = round(Ncomm * parameters$probGrVar)
groupVar = rep(0, Ncomm)
groupVar[sample(Ncomm, groupsWithEffect)] = 1

#groupVar =  rbinom(parameters$Ncomm, 1, parameters$probGrVar)
#noGroupEff = groupVar==0

indiVar= rbinom(parameters$size, 1, parameters$probIndiVar)
#noIndiVar = indiVar==0
#missIndiVar = rbinom(parameters$size,1, parameters$IndiMissRate) & !noIndiVar
#indiVar[missIndiVar]=1


result$groupVar = groupVar[result$group]
result$indiVar = indiVar


regionRandomEff = rnorm(Ncomm, 0, parameters$regionsd)


probsLogit = log(parameters$mu)+regionRandomEff +
             log(parameters$DiffEffSizeGroup)*groupVar+
             log(parameters$DiffEffSizeIndi)*indiVar+
             log(parameters$DiffEffSizeInter)*groupVar*indiVar
             
probs = exp(probsLogit) / (1+exp(probsLogit))             


smoke = rbinom(length(probs), 1, probs)
# misclassfy smoke
#noSmoke = smoke==0
#missSmoke = rbinom(parameters$size,1, parameters$SmokeMissRate) & !noSmoke
#smoke[missSmoke]=1


result$smoke=smoke

return(result)


}

#############################################
#generate new dataset and calculate SE #
#############################################                                                             
SEvalue=function(Ngroups, parameters, varying,nullmodel,
  Nmethods) {
  
estparaR= c("intercept","groupvar","indivar","interaction","interceptRan")
estparaC = c("est","sd","LoweCI", "UpperCI")

if (nullmodel==T)  {
varying=c(varying,list("DiffEffSizeGroup"=1, "DiffEffSizeIndi"=1,"DiffEffSizeInter"=1))
estparaR = c("intercept", "interceptRan")
}
EstCI=array(NA, c(length(estparaR),length(estparaC),length(Nmethods),length(Ngroups)))
dimnames(EstCI)=list(estparaR,estparaC,Nmethods,paste("Ngroups", "=" , Ngroups, sep=""))



for (Dgroup in 1:length(Ngroups)) {

subdata = simCohort(parameters, Ncomm=Ngroups[Dgroup], varying)

if(any(Nmethods == "lmer")) {
detach("package:nlme")
#res3=NULL
if (nullmodel==T) {res3 = try(lmer(smoke ~ (1|group), data=subdata, family='binomial'))
}else{                                                
res3 = try(lmer(smoke ~ groupVar + indiVar + groupVar*indiVar+(1|group), data=subdata, family='binomial'))
}
EstCI[,,"lmer",Dgroup]=matrix(rep("NA",20),nrow=length(estparaR), ncol=length(estparaC))
if(is.null(names(res3))) {
allEstlmer = c(fixef(res3), "sigma"=sqrt(as.numeric(VarCorr(res3)$group))) 
allSElmer = c(sqrt(diag(vcov(res3))), NA)

EstCI[,,"lmer",Dgroup]=cbind("est"=allEstlmer, 
                             "sd"=allSElmer, 
                             "LowCI"=allEstlmer-1.96* allSElmer, 
                             "UpperCI"=allEstlmer+1.96*allSElmer) 

} 

}

library(nlme)
library(R2WinBUGS)
library(glmmBUGS)
forBugs = NULL

if (any(Nmethods=="mcmc"|Nmethods=="pql")) {

if (nullmodel == T) {
#try(
forBugs = glmmBUGS(smoke ~ 1 , effects="group",  family='bernoulli', data=subdata)
modelFile = "modelNull.bug"
#)
}else{
#try(
forBugs = glmmBUGS(smoke ~ groupVar + indiVar + groupVar*indiVar , effects="group", family='bernoulli', data=subdata)
modelFile = "modelFull.bug"    
forBugs$startingValues$betaobservations[1:2]= 0
forBugs$startingValues$betagroup= 0
                                                                           
#)
}
res2 = forBugs$pql

EstCI[,,"pql",Dgroup]=matrix(rep(NA,length(estparaR)*length(estparaC)),nrow=length(estparaR), ncol=length(estparaC))
if(!is.null(names(res2))) {
allEstpql = as.numeric(c(summary(res2)$tTable[,1],"sigma"=sqrt(as.matrix(res2$modelStruct$reStruct$group)) * res2$sigma))

allSEpql = c(as.numeric(summary(res2)$tTable[,2]),NA)


EstCI[1:length(allEstpql),,"pql",Dgroup]=cbind("est"=allEstpql, 
                "sd"=allSEpql, 
                "LowCI"=c(allEstpql-1.96* allSEpql), 
                "UpperCI"=c(allEstpql+1.96*allSEpql)) 
}
}

############ winbugs here ###############

if(any(Nmethods == "mcmc"  )) {
forBugs$startingValues$vars$group  = 0.3
forBugs$startingValues$Rgroup = rnorm(length(forBugs$startingValues$Rgroup), 
  forBugs$startingValues$intercept, 0.1)

assign("testStuff", forBugs, pos=1)
                     
assign("startingValues", forBugs$startingValues, pos=1)




source("getInits.R")
fromBugs = bugs(forBugs$ragged, getInits, 
  parameters.to.save = names(getInits()),
  model.file=modelFile, n.chain=3, n.iter=3100, 
  n.burnin=100, n.thin=10, 
	working.directory="/home/hedy/cohortsamplesize/smoke", 
	save.history=F) 
print(1)
bugsResult = summaryChain(restoreParams(fromBugs, forBugs$ragged))
print(2)
bugsResult =  rbind(bugsResult$scalars, bugsResult$betas)
print(3)
# merge results in
if (nullmodel==T) {
EstCI[,,"mcmc",Dgroup] = 
  bugsResult[
    c("intercept","SDgroup"),
    c("mean","sd","2.5%","97.5%")
  ]
}else{
EstCI[,,"mcmc",Dgroup] = 
  bugsResult[
    c("intercept","group","indiVar","groupVar_indiVar","SDgroup"),
    c("mean","sd","2.5%","97.5%")
  ]
}
}


############ change to binomial ##########

Count=binToBinom(subdata$smoke, subdata[,c("groupVar","indiVar","group")])
subdata$ones = 1
CountGroup =binToBinom(subdata$smoke, subdata[,c("group", "ones")])


if(length(grep("^gee", Nmethods))) {
############  GEE starts here ############
if (nullmodel == T) {
  
resgee = try(gee(cbind(y,N-y)~ 1, id=group, data=CountGroup, family=binomial, corstr="exchangeable")
)
 
}else{
resgee = try(gee(cbind(y,N-y)~ groupVar + indiVar + groupVar*indiVar, nid=group, data=Count, family=binomial, corstr="exchangeable")
)}
if(!is.null(names(resgee))) {
geeEst = resgee$coefficients
geeSErobust = sqrt(diag(resgee$robust.variance))
geeSEnaive = sqrt(diag(resgee$naive.variance))
EstCI[,,"geeRobust",Dgroup]=rbind(cbind("est"=geeEst, 
                "sd"=geeSErobust, 
                "LowCI"=c(geeEst-1.96* geeSErobust), 
                "UpperCI"=c(geeEst+1.96*geeSErobust)), rep(NA,dim(EstCI)[2])) 


EstCI[,,"geeNaive",Dgroup]=rbind(cbind("est"=geeEst, 
                "sd"=geeSEnaive, 
                "LowCI"=c(geeEst-1.96* geeSEnaive), 
                "UpperCI"=c(geeEst+1.96*geeSEnaive)),rep(NA,dim(EstCI)[2])) 
}else{
EstCI[,,"geeRobust",Dgroup]=EstCI[,,"geeNaive",Dgroup]=NA
}
}

############## INLA starts from here #######
if(any(Nmethods == "inla")){
write.table(Count, file="count.txt")
if (nullmodel == T)    {
formula=y~f(group, model="iid", param=c(.001, 0.001))
} else {
Count$interaction = Count$groupVar * Count$indiVar
formula=y~groupVar+indiVar+interaction+f(group, model="iid", param=c(.001, 0.001) )
}
resInla=NULL
resInla = try(inla(formula,data=Count,family="binomial",Ntrials=N,
	quantiles = c(0.025,0.975), control.inla=list(dz=0.25)))

if(class(resInla)!="try-error") {
	resSEran =  getInlaSDsummary(resInla)


# fill in fixed effects
	qq=length(rownames(resInla$summary.fixed)) 

	EstCI[1:qq,,"inla", Dgroup] =
		resInla$summary.fixed[,-5]
# fill in sd
	EstCI["interceptRan",c(names(resSEran)), "inla", Dgroup] = resSEran
} # end if !is.null

}
}
result=list()
result$EstCI = EstCI
result$subdata = subdata
return(result)
}

###### add simulation time #########
allEstCI = function(Ngroups,parameters, varying,nullmodel,simulationTime, Nmethods = c("lmer","pql","geeNaive", "geeRobust","inla")) {
 
 EEEstCI = SEvalue(Ngroups, parameters, varying, nullmodel,Nmethods)
 EstCI = EEEstCI$EstCI
 subdata= EEEstCI$subdata
 allEstCI = array(NA, c(dim(EstCI), simulationTime))
 dimnames(allEstCI) = c(dimnames(EstCI), list(as.character(1:simulationTime)))
 allEstCI[,,,,1]=EstCI
 
 #write.csv(subdata, file=paste(varying, nullmodel, "1",".txt"))
 for (i in 2:simulationTime){
   EEEstCI = SEvalue(Ngroups, parameters, varying, nullmodel, Nmethods) 
   EstCI = EEEstCI$EstCI
   subdata= EEEstCI$subdata
   allEstCI[,,,,i] = EstCI
  
# write.csv(subdata, file=paste(varying, nullmodel, i,".txt"))
 }

 if (is.character(allEstCI)) 
 allEstCI= type.convert(allEstCI)
 
 return(allEstCI)
}


################# calculate coverage and bias####################
CoverBias = function( allEstCI, parameters, varying, nullmodel, simulationTime) {
# calculate the non-coverage and bias
Nmethods = dimnames(allEstCI)[[3]]
Ngroups = as.numeric(gsub("^Ngroups=","",dimnames(allEstCI)[[4]]))
simulationTime=dim(allEstCI)[5]

result=NULL
if (nullmodel==T)  {
varying=c(varying,list("DiffEffSizeGroup"=1, "DiffEffSizeIndi"=1,"DiffEffSizeInter"=1))
}

for (Dvarying in names(varying))
    parameters[[Dvarying]] = varying[[Dvarying]][1]

if (nullmodel == T) {
TruePara = c("intercept"=log(parameters$mu), "interceptRan"=parameters$regionsd)
} else{
TruePara = c("intercept"=log(parameters$mu), "groupvar"=log(parameters$DiffEffSizeGroup), "indivar"=log(parameters$DiffEffSizeIndi),
"interaction"=log(parameters$DiffEffSizeInter),"interceptRan"=parameters$regionsd)
}

nNAintercept = apply(is.na(allEstCI[1,1,,,]),c(1,2),sum)
            
Noncoverage =Covera = biasE = array(NA, c(dim(allEstCI)[1], dim(allEstCI)[3:4]))
dimnames(Noncoverage) =dimnames(biasE)= dimnames(Covera)= c(list(dimnames(allEstCI)[[1]]), list(dimnames(allEstCI)[[3]]), list(dimnames(allEstCI)[[4]]))
for (Dest in 1:length(names(TruePara))) {
     for (Dmethod in 1: dim(allEstCI)[3]) {
         for (DNgroup in 1:length(Ngroups)){
            Ntot = simulationTime - nNAintercept[Dmethod,DNgroup]

            AllCI= allEstCI[names(TruePara)[Dest], c("LoweCI","UpperCI"),Dmethod,DNgroup,]
            Noncoverage[Dest,Dmethod,DNgroup]=sum(!(AllCI[1,]< TruePara[Dest] & TruePara[Dest]< AllCI[2,]), na.rm=T )
            Covera[Dest,Dmethod,DNgroup]= 1 - (Noncoverage[Dest,Dmethod,DNgroup]/Ntot)
            AllEstimate= allEstCI[names(TruePara)[Dest], "est",Dmethod,DNgroup,]
            biasE[Dest,Dmethod,DNgroup]=(mean(AllEstimate, na.rm=T) - TruePara[Dest])/TruePara[Dest]*100
}
}
}

#Covera["interceptRan", c("lmer","pql","geeNaive","geeRobust"),]=NaN

result$nNAintercept=nNAintercept
result$Covera=Covera
result$biasE = biasE
return(result)                                  
}


##############################design effect #######################

fdeff= function (parameters, Ngroups, varying){
for (Dvarying in names(varying))
    parameters[[Dvarying]] = varying[[Dvarying]][1]

seDEFF=matrix(NA, nrow=4, ncol=length(Ngroups), dimnames=list(c("DEFFsim","DEFFapprox","SEsim","SEapprox"),c(as.character(Ngroups))))

meanOnLogitScale = log(parameters$mu)

for (Dgroup in 1:length(Ngroups)) {
DEFF = getDEFF(meanOnLogitScale, sigma=parameters$regionsd, 
  Nsim = 5000,Nsubj=parameters$size/Ngroups[Dgroup] )

#logitSE = 1/sqrt(parameters$size/DEFF*parameters$mu/((1+parameters$mu)^2))
# I think it should be this instead?
meanNaturalScale = attributes(DEFF)$moments[1]
logitSE = 1/sqrt( (parameters$size/DEFF) * meanNaturalScale*(1-meanNaturalScale))


seDEFF[,Dgroup]=c(DEFF, logitSE)
}
return(seDEFF)
}

#############design effect calculation from here##########################
getDEFF = function(
    logitMean, sigma, 
    Nsim = 200000,
    Nsubj ) {

  # simulate some normals to find the moments of the logit-normal  
  theSim = rnorm(Nsim, logitMean, sigma)
  # transform to logit-normals
  theSim = exp(theSim) / (1 + exp(theSim))
  # compute the moments
  moment1 = mean(theSim)
  moment1sq = moment1^2
  moment2 = mean(theSim^2)
  
  # the ICC
 # c(moments=(moment2 - moment1sq) / (moment1 - moment1sq),
 #   corr = cor(rbinom(Nsim, 1, theSim), rbinom(Nsim, 1, theSim)) )
  
 ICC = c(sim=(moment2 - moment1sq) / (moment1 - moment1sq),
            calc = sigma^2 / (sigma^2 + pi^2/3) )

 result=1+(Nsubj-1)*ICC
  names(result)=c("simulation","approx")
  
  attributes(result)$moments = c(moment1, moment2)
  
#if(!is.null(Nsubj)) result = 1+(n-1)*ICC

result
 }

getInlaSDsummary = function(resInla) {
	sdmargin = resInla$marginals.hyperpar[[1]]
	colnames(sdmargin)=c("x", "y")
	sdmargin = cbind(sdmargin, sdX=NA, diffx = NA)
	sdmargin[,"sdX"] = sqrt(1/sdmargin[,"x"])
	sdmargin[,"diffx"] = c(diff(sdmargin[,"x"]),0) 	
	result = c(est=sum(sdmargin[,"y"]*sdmargin[,"sdX"]*sdmargin[,"diffx"]))
	result["sd"] =
	 sqrt(sum(sdmargin[,"y"]*sdmargin[,"sdX"]^2*sdmargin[,"diffx"]) - result["est"]^2)
	
	cumulativeProb =cumsum(sdmargin[,"y"] * sdmargin[,"diffx"])			
	cumulativeProb = cumulativeProb/max(cumulativeProb)
	
	result[c("LoweCI","UpperCI")] = 
		sort(sdmargin[c(
			order(abs(cumulativeProb-0.025), decreasing=F)[1],
			order(abs(cumulativeProb-0.975), decreasing=F)[1]),
			     "sdX"])
	result
}
             
SEgraph = function (allEstCI, CovBia, nullmodel, graphName) {
seIntercept = apply(allEstCI["intercept","sd",,,],c(1,2),function(x) mean(x,na.rm=T))
seIntercept = t(rbind(seIntercept, CovBia$seDEFF[c("SEsim","SEapprox"),]))
if (nullmodel==F) {
seGroupVar = t(apply(allEstCI["groupvar","sd",,,],c(1,2),function(x) mean(x,na.rm=T)))
seIndiVar = t(apply (allEstCI["indivar","sd",,,], c(1,2),function(x) mean(x, na.rm=T)))
seInteraction = t(apply(allEstCI["interaction", "sd",,,], c(1,2),function(x) mean(x, na.rm=T)))  
}


theLables=gsub("[[:alnum:]|:]+=","",rownames(seIntercept))
theLables=gsub(",","",theLables)
theLables=as.numeric(theLables)
pdf(graphName)
Scolour=c("black","blue","green", "grey", "red","orange", "purple")
Slty = rep(1,7)
par(mfrow=c(2,2))

matplot(theLables, seIntercept, col=Scolour, 
 xlim=c(1,max(theLables)),ylim=c(min(seIntercept,na.rm=T), max(seIntercept,na.rm=T)),
 lwd=c(1,1,2,2,1),lty=Slty, 
 xlab="The Number of Groups", ylab ="Standard Error of Intercept", pch=1, type="o")

if (nullmodel==F){
matplot(theLables, seGroupVar, col=Scolour, 
 xlim=c(1,max(theLables)),ylim=c(min(seGroupVar,na.rm=T),max(seGroupVar,na.rm=T)), 
 lwd=c(1,1,2,2,1),lty=Slty, 
 xlab="The Number of Groups", ylab ="Standard Error of Group Level Variable", pch=1, type="o")

matplot(theLables, seIndiVar, col=Scolour, 
 xlim=c(1,max(theLables)),ylim=c(min(seIndiVar, na.rm=T),max(seIndiVar,na.rm=T)),
 lwd=c(1,1,2,2,1),lty=Slty, 
 xlab="The Number of Groups", ylab ="Standard Error of Individual Level Variable", pch=1, type="o")

matplot(theLables, seInteraction, col=Scolour, 
 xlim=c(1,max(theLables)),ylim=c(min(seInteraction,na.rm=T),max(seInteraction,na.rm=T)), 
 lwd=c(1,1,2,2,1),lty=Slty, 
 xlab="The Number of Groups", ylab ="Standard Error of Group-Individual Level Interactions", pch=1, type="o")
}
legend("topright",colnames(seIntercept), col=Scolour, lty=Slty)

dev.off()
}

CoverBiasgraph = function (allEstCI, Covera, biasE, nullmodel, CgraphName, BgraphName) {

Cintercept = t(Covera["intercept",,])
Bintercept = t(biasE["intercept",,])

BinterceptRan = t(biasE["interceptRan",c("lmer","pql","inla"),])
if (nullmodel==F) {
Cgroupvar = t(Covera["groupvar",,])
Bgroupvar = t(biasE["groupvar",,])
Cindivar = t(Covera["indivar",,])
Bindivar = t(biasE["indivar",,])
Cinteraction = t(Covera["interaction",,])
Binteraction = t(biasE["interaction",,])
}

theLables=gsub("[[:alnum:]|:]+=","",rownames(Cintercept))
theLables=gsub(",","",theLables)
theLables=as.numeric(theLables)
Scolour=c("black","blue","green", "grey", "red","orange", "purple")
Slty = rep(1,7)

pdf(CgraphName)

par(mfrow=c(2,2))
matplot(theLables, Cintercept, col=Scolour, 
 xlim=c(1,max(theLables)),ylim=c(min(Cintercept,na.rm=T), max(Cintercept,na.rm=T)),
 lwd=c(1,1,2,2,1),lty=Slty, 
 xlab="The Number of Groups", ylab ="Coverage of Intercept", pch=1, type="o")

if (nullmodel==F){
matplot(theLables, Cgroupvar, col=Scolour, 
 xlim=c(1,max(theLables)),ylim=c(min(Cgroupvar,na.rm=T), max(Cgroupvar,na.rm=T)), 
 lwd=c(1,1,2,2,1),lty=Slty, 
 xlab="The Number of Groups", ylab ="Coverage of Group Level Variable", pch=1, type="o")

matplot(theLables, Cindivar, col=Scolour, 
 xlim=c(1,max(theLables)),ylim=c(min(Cindivar,na.rm=T), max(Cindivar,na.rm=T)),
 lwd=c(1,1,2,2,1),lty=Slty, 
 xlab="The Number of Groups", ylab ="Coverage of Individual Level Variable", pch=1, type="o")

matplot(theLables, Cinteraction, col=Scolour, 
 xlim=c(1,max(theLables)),ylim=c(min(Cinteraction,na.rm=T), max(Cinteraction,na.rm=T)),
 lwd=c(1,1,2,2,1),lty=Slty, 
 xlab="The Number of Groups", ylab ="Coverage of Group-Individual Level Interactions", pch=1, type="o")
}
legend("topright",colnames(Cintercept), col=Scolour, lty=Slty)

dev.off()


pdf(BgraphName)

par(mfrow=c(3,2))

matplot(theLables, Bintercept, col=Scolour, 
 xlim=c(1,max(theLables)),ylim=c(min(Bintercept,na.rm=T), max(Bintercept,na.rm=T)),
 lwd=c(1,1,2,2,1),lty=Slty, 
 xlab="The Number of Groups", ylab ="Bias of Intercept", pch=1, type="o")
legend("topright",colnames(Bintercept), col=Scolour, lty=Slty)

matplot(theLables, BinterceptRan, col=Scolour, 
 xlim=c(1,max(theLables)),ylim=c(min(BinterceptRan,na.rm=T), max(BinterceptRan,na.rm=T)),
 lwd=c(1,1,2,2,1),lty=Slty, 
 xlab="The Number of Groups", ylab ="Bias of Random Effect", pch=1, type="o")
legend("topright",colnames(BinterceptRan), col=Scolour, lty=Slty)

if (nullmodel==F){
matplot(theLables, Bgroupvar, col=Scolour, 
 xlim=c(1,max(theLables)),ylim=c(min(Bgroupvar,na.rm=T), max(Bgroupvar, na.rm=T)),
 lwd=c(1,1,2,2,1),lty=Slty, 
 xlab="The Number of Groups", ylab ="Bias of Group Level Variable", pch=1, type="o")
legend("topright",colnames(Bgroupvar), col=Scolour, lty=Slty)

matplot(theLables, Bindivar, col=Scolour, 
 xlim=c(1,max(theLables)),ylim=c(min(Bindivar,na.rm=T), max(Bindivar, na.rm=T)),
 lwd=c(1,1,2,2,1),lty=Slty, 
 xlab="The Number of Groups", ylab ="Bias of Individual Level Variable", pch=1, type="o")
legend("topright",colnames(Bindivar), col=Scolour, lty=Slty)

matplot(theLables, Binteraction, col=Scolour, 
 xlim=c(1,max(theLables)),ylim=c(min(Binteraction, na.rm=T), max(Binteraction, na.rm=T)), 
 lwd=c(1,1,2,2,1),lty=Slty, 
 xlab="The Number of Groups", ylab ="Bias of Group-Individual Level Interactions", pch=1, type="o")
legend("topright",colnames(Binteraction), col=Scolour, lty=Slty)

}
dev.off()
}


